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Abstract 

Two-component Fermi gases with tunable repulsive or attractive interactions inside quasi-one- 
dimensional (Q1D) harmonic wells may soon become the cleanest laboratory realizations of strongly 
correlated Luttiger and Luther-Emery liquids under confinement. We present a microscopic Kohn- 
Sham density- functional theory of these systems, with specific attention to a gas on the approach to 
a confinement-induced Feshbach resonance. The theory employs the one-dimensional Gaudin-Yang 
model as the reference system and transfers the appropriate Q1D ground-state correlations to the 
confined inhomogeneous gas via a suitable local-density approximation to the exchange and cor- 
relation energy functional. Quantitative understanding of the role of the interactions in the bulk 
shell structure of the axial density profile is thereby achieved. While repulsive inter component 
interactions depress the amplitude of the shell structure of the noninter acting gas, attractive in- 
teractions stabilize atomic-density waves through spin pairing. These should be clearly observable 
in atomic clouds containing of the order of up to a hundred atoms. 



PACS numbers: 03.75.Ss,71.15.Mb,71.10.Pm 



I. INTRODUCTION 



Ultracold atomic gases, which are highly tunable and ideally clean, are attracting a great 
deal of interdisciplinary interest. In particular their study may help us understand a number 
of phenomena that have been predicted in solid-state and condensed- matter physics jl[ . Sev- 
eral effects, known in these subfields of physics for decades, have already been observed and 
quantitatively analyzed in ultracold atomic gases. Three beautiful examples are the Bloch 
oscillations under an applied force in a one- dimensional (ID) optical lattice 0, the forma- 
tion of highly ordered Abrikosov lattices of vortices in rapidly rotating harmonic traps 2| , 
and the superffuid-to-Mott insulator transition of a condensate in a 3D optical lattice |4|. 

Cold atoms have also been successfully trapped in low-dimensional geometries . Typical 
ID quantum phenomena have already been observed in both Bose and Fermi gases. For 
instance, in the work of Paredes et al. and of Kinoshita et al. ||| a 87 Rb gas has been 
used to realize experimentally a Tonks- Girardeau system |f|. The more recent preparation of 
two-component Fermi gases in a quasi-l-D (QID) geometry [?J provides a unique possibility 
to experimentally study phenomena that were predicted a long time ago for electrons in a 
ID solid-state environment, such as spin-charge separation in Luttinger liquids [H, 9| and 
charge- density waves in Luther- Emery liquids [10j . The experiment by Moritz et al. 71 also 
offers the opportunity of testing a Q1D integrable model of the BCS-BEC crossover [ll. 12l|. 
which is based on the idea of a confinement- induced resonance (CIR) |l3| |. 

In 11, l3| the gas has been assumed to be translationally invariant along the axial 
direction, and thus the authors have been able to provide an analytical description of the 
crossover by employing the exact Bethe-Ansatz (BA) solution of the homogenous Gaudin- 
Yang model (HGYM) 0] and of the homogeneous Lieb-Liniger model The present 

work focuses instead on inhomogeneous Q1D Fermi gases inside highly elongated harmonic 
traps and treats the axial confinement by means of the Hohenb er g-K ohn- Sham density- 
functional theory (DFT) HQ. With a few exceptions [l7L llH. Hoi 120^ . most applications 
of DFT have used so far as the underlying reference fluid the homogeneous electron gas, 
which is a normal Fermi liquid over a wide range of density. In our present study we use the 
HGYM as the reference fluid, in order to transfer to the inhomogeneous gas the Luttinger 
and Luther-Emery ID correlations. 

It is appropriate at this point to refer to related theoretical studies dealing with Q1D 
inhomogeneous Fermi eases In Ref. 21] a bosonization technique has 



been used to calculate analytically the density profile, the momentum distribution, and sev- 
eral correlation functions of two-component Fermi gases with inclusion of intercomponent 
forward-scattering processes. In Refs. [13, 13, 24[ the Thomas-Fermi approximation (see 



the discussion in Sect. Ill Al below) and the so-called inhomogeneous Tomonaga-Luttinger 
liquid model have been used to calculate the density profile of a large system and to discuss 
spin-charge separation in two-component Fermi gases. In the present work we perform mi- 
croscopic calculations of the ground-state (GS) density profile of systems with arbitrary size, 
without having to assume neither peculiar intercomponent interactions (as in the bosoniza- 
tion scheme of Ref. j2l|) nor very large atom numbers (as in Refs. |22l 0, 0|). We give a 
fully quantitative study of how exchange and correlations modify the bulk shell structure of 
the axial density profile. In particular we show that for sufficiently strong attractive interac- 
tions, experimentally detectable atomic-density waves (ADWs) are formed by spin pairing 
along the axial direction, which should be clearly observable in systems with a relatively low 
number of atoms (iVf < 100). These the oversimplified Thomas- Fermi treatments cannot 
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predict. 

The contents of the paper are briefly as follows. In Sect. |H]we introduce the Hamilto- 
nian that we have used to describe the system of present physical interest, summarize the 
properties of the model in the absence of external potentials, and describe the self-consistent 
DFT scheme that we have used to deal with the inhomogeneity. In Sect. IHII we report and 
discuss our main numerical results and finally in Sect. IIVI we draw our main conclusions. 
An Appendix contains the exact solution of the inhomogeneous model for two atoms only, 
which is used in the main text for a test of the local density approximation in the extreme 
limit of low particle numbers. 



II. THEORETICAL APPROACH 

We consider a two-component Fermi gas with N{ atoms confined inside a strongly elon- 
gated harmonic trap. The two species of fermionic atoms are assumed to have the same 
mass m and different pseudospin a, a =] or [. The trapping potential is axially symmetric 
and is characterized by angular frequencies u± and uin in the radial and the longitudinal 
directions, with uuu <C u±. Correspondingly we introduce the harmonic-oscillator lengths 
a± = y/hj (muj±) and a\\ = JKj (mu>\\). 

The gas is dynamically ID if the anisotropy parameter of the trap is much smaller than 
the inverse atom number, u\i/luj_ <C Nf 1 . It can thus be described by the Hamiltonian |24| 

h 2 Ni d 2 Nl Nl 

i=l 1 i=l j=l 

neglecting intracomponent p-wave interactions. Here, 

= 2h 2 a 3D (B) 1 
1D ma 2 ^ 1 — Aa 3U (B)/a± 



is the effective ID Olshanii coupling parameter [13j, with A = |C(l/2)|/v / 2 ~ 1.0326 and 

being the Hurwitz zeta function, and V ext = Y^i=iVext( z i) = ( mu, f\/^) ^2^=1 z f * s ^ ne 
external static potential associated with the axial confinement. The 3D scattering length 
a 3D can be tuned by means of a magnetic field B and has the resonant structure a 3D (B) = 
ab g [l — SB/(B — Bp)], Bp being the position of a Feshbach resonance, 5B its width, and a^g 
the so-called background scattering length 0]. 

Choosing the harmonic-oscillator length an as unit of length and the harmonic-oscillator 
quantum fkov, as unit of energy, the Hamiltonian (0) can be shown to be governed by the 
dimensionless coupling parameter 

A = nf- (3) 

The parameter A diverges at the CIR, i.e. when the external magnetic field takes the value 
B* = Bp — 5B(a* D /ab g — with a* D = a±/A. The coupling parameter is negative for 
B* < B < Bp + 5B and positive everywhere else. At the 3D Feshbach resonance, i.e. when 
B = Bp, the coupling parameter has the finite value A F = — 2a||/(^4aj_). In Fig. [T] we show 
the dependence of A on the magnetic field B. 
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For B > B* and V ext = the Hamiltonian (0) reduces to the HGYM, which can be solved 
exactly by means of the BA technique for both repulsive (g m > 0) and attractive (g 1D < 0) 
interactions [14J. In the thermodynamic limit (iVf, L — > oo, L being the system size) and for 
a pseudosp in-compensated system (iVj = N±), the properties of the HGYM are determined 
by the linear density n = Nf/L and by the effective coupling g 1D . These can be conveniently 
combined into a single dimensionless parameter 7 = mg lu /(h 2 n). 

The energy per atom can be written in terms of the "momentum distribution" p(k) as 

. e h 2nu f +Q dk h 2 k 2 ns 

£ GS {n, g lu ) = — + / — - — p(k) , 4 

2 n J _q 2tt 2m 

where £b = 0, v — 1 for g 1D > and £ b = —mg 2 D / (4^, 2 ), v = 2 for g 1D < 0. The function p{k) 
can be calculated by solving the Gaudin-Yang BA integral equation jlj 



P(k) = f + -[ +Q p-JC (2(k - q )/{ in )) p{q) , (5) 
2tt 7n J_ Q 2tt 



where Q is determined by the normalization condition 



+Q n 

p{k)dk = - (6) 

n V 



and the kernel JC(x) is given by 



30 sech Try/2 ^ n 
AC(x)= \ J-^ \ l + ^ + v) 2 ] . (7) 

TT^Vi for ^ <0 

For (jf 1D > the HGYM describes a Luttinger liquid 0, 0|, while for g m < it describes a 
Luther- Emery liquid [lo| . 

Before proceeding to discuss the properties of the inhomogeneous model under confine- 
ment, we should stress that Eq. (0))-© describe the homogeneous limit of the model only 
for B > B*, i.e. before the CIR. After the CIR, as discussed in Refs. [lli Il2|. the fermion 
pairs become unbreakable spin-singlet dimers, behaving like bosons with mass 2m and den- 
sity n/2. Thus the appropriate homogeneous limit for B < B* is the Lieb-Liniger gas of 



ap 

interacting bosons 15| , and one should resort in treating inhomogeneity to a DFT approach 



such as that proposed by Griffin |17( (see also the work of Oliveira et al. |17|). 



A. Density- functional theory of QID gas in the Kohn-Sham scheme 

In the presence of a longitudinal external potential the Hamiltonian 7i in Eq. (jlj cannot 
be diagonalized exactly We calculate the GS properties of TC for B > B* by resorting to 
the fermionic DFT scheme 0, 0| • 

Within the Kohn-Sham version of DFT the GS density, n GS (z) = (GS | J2 t 8(z — Zj)|GS), 
can be calculated by solving self-consistently the Kohn-Sham-Schrodinger (KSS) equa- 
tions llfil . 

ip a (z) = e a <f a (z) (8) 



2m dz 2 



+ Vks(,z; [n GS (z)}) 
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with VksO; [n G s(z)}) = v u (z; [n GS (z)]) + v xc (z; [n GS (z)}) + V cxt (z), together with the closure 

n GS (z) = ^2 Ta l^«( z )| 2 ■ (9) 

o,occ. 

Here the sum runs over the occupied orbitals and the degeneracy factors T a satisfy the sum 
rule ^ Q T a = Nf . The first term in the effective Kohn-Sham potential Vks is the Hartree 
term vu = gmn GS (z), while the second term is the exchange-correlation (xc) potential defined 
as the functional derivative of the xc energy £ xc [n(z)] evaluated at the GS density profile, 
v xc = SS xc [n(z)]/Sn(z)\ GS . The total GS energy of the system is given by 

/+oo 
dzv xc (z; [n GS (z)])n GS (z) 
-oo 

i 

n 2 GS (z)dz + £ xc [n GS (z)] . (10) 
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ID / ^2 

2 



oo 



Equations (jSJ) and JHJ) provide a formally exact scheme to calculate n GS (z) and S GS , but 
S xc and v xc need to be approximated. The local-density approximation (LDA) has been 
shown toprovide an excellent description of the GS properties of a variety of inhomogeneous 
systems jil[l^. In the following we employ a BA-based LDA (BALDA) functional |l9l Eol] 
for the xc potential, 

vir\^ K s (z)}) = v^(n,g m )\ n _> nGs[z) . (11) 

Here the xc potential of the HGYM is defined by 

d 

w x° m ( n > 9m) = [ne GS {n, g m ) - nK(n)} - ng 1D , (12) 

K,(n) = 7T 2 h 2 n 2 /(24m) being the kinetic energy of the noninteracting gas per atom. 

Before discussing specific calculations of the xc potential of the HGYM, several important 
remarks are in order at this point: 

(i) In the limit A = the KSS equations correctly yield the GS density profile of a 
noninteracting paramagnetic Fermi gas, 

n GS (z)\ x=0 = — ^exp(-z /o||) ^ 2 » n , , ( 13 ) 

II * n=0 

given in terms of the Hermite polynomials H n (x) of degree < n < N{/2 — 1. This density 
profile exhibits a shell structure characterized by iVf/2 oscillations, whose origin lies in the 
fermionic statistical correlations: the occupation probability V(n) = ^2 a (cn a Cn,a) of the 
ID harmonic-oscillator states is unity for < n < Nf/2 — 1 and zero for n > Nf/2. The 
existence of this sharp "Fermi edge" is ultimately responsible for the bulk shell structure, 
which is analogous to the Friedel oscillations originating in a normal Fermi liquid from the 
sharply defined Fermi surface 0. The occupation probabilities and the shell structure are 
expected to be strongly affected by many-body xc effects (see Sect. IIHI below). 
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(ii) In the limit A = +00 the GS density profile should become that of a fully spin- 
polarized noninteracting Fermi gas, 

n GS (z)\ x=+oo = -i-exp(-z 2 /4) Hn 2*ri^ > (14) 

II V n=0 

exhibiting a shell structure characterized by Nf oscillations [Hf. This asymptotic property 
can be checked explicitly for Nf — 2 (see the Appendix) and originates from the fact that an 
infinitely strong repulsion between antiparalleL-pseudospin atoms in ID acts like the Pauli 
principle between parallel-pseudospin atoms |22|, |24| . The present formalism does not apply 
to such a strong coupling regime (note also that for B > B* the value of A is bounded from 
above). 

(iii) The main difference between the present BALDA scheme and the Thomas-Fermi 
approach is that in the latter (22! . EiL E^ the LDA is also used to approximate the nonin- 
teracting kinetic energy functional T s [n GS (z)], which is written as 

/+00 
[nK(n)] n ^ nas{z) dz . (15) 
-00 

The (Hohenberg-Kohn) Thomas-Fermi equation reads 

+ V KS (z; [n GS (z))) = constant , (16) 

n-m G s(*) 

the constant being fixed by normalization. The Thomas-Fermi profile misses the shell struc- 
ture as well as atom tunnel beyond the Thomas- Fermi radius Z^f- In our approach, instead, 
Tslrics^z)} is treated exactly through the Kohn-Sham mapping 

h 2 f +OQ d 2 
T s [n GS (z)] = -— J ^ ^ a (z)— Va (z)dz , (17) 

the single-particle orbitals ip a (z) = f a (z] [n GS (z)}) being unique functionals of the GS den- 
sity UM- 



d[nK(n)] 
dn 



B. The exchange-correlation potential 



In what follows we propose two different ways to calculate the xc potential of the HGYM. 



1. BALDA/1 

The potential v^° m (n, g m ) can be calculated by applying its definition (|T2j) directly to 
Eqs. (Jlj)-©- It is easy to show that the xc potential of the HGYM is exactly given by the 
following equation, 



£ b „ &Q 2 dQ _ f +Q dk h 2 k 2 dp{k) 



■2 



.-(„, gm) = f + 2^p(Q) S + 2 m J £ - - ngm , (18) 
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where d n Q and d n p satisfy the coupled BA equations 

= ^l^(^-Q)/(in))+}C(2(k + Q)/( 1 n))]^ 
v 



and 



+ TnL d £ K(2(k - k ' )l(in))e -^r < 19 > 



« /^N^Q f +Q dp(k) „ 1 



An accurate numerical solution of these coupled BA equations leads to the exact xc potential 
of the HGYM. 

The results for n GS (z) that are obtained with v^° m (n, g 1D ) determined according to this 
route will be termed with the acronym BALDA/1. 



2. BALDA/2 

As an alternative v^° m (n, g 1D ) can also be calculated from accurate analytical parametriza- 
tions of the GS energy of the HGYM, which incorporate exactly known limiting behaviors 
both at weak and strong coupling. This route will reduce the numerical effort and affords a 
test of the sensitivity of the results to the details of the implementation of the theory. 

Let us introduce the Fermi wave number kp = ixn/2 and the Fermi energy Bp = 
h 2 kp/(2m). For repulsive interactions we find that the GS energy of the HGYM in units of 
the Fermi energy, e = e GS /ep, can be very accurately parametrized by the simple formula 

, . 4x 2 /3 + a v x + b„ 
e(x>0) = — { p p , (21) 

where x = a p = 5.780126, b p = -(8/9) In 2 + twJA, c p = (8/vr)ln2 + 3a p /4, and 

dp = 3b p . Equation (|2*T|) embodies the exact behaviors |24j 

. 1/3 + x/tt + ... fors^0+ / N 

e(x)={ . (22) 

4/3 - 321n2/(37rx) + ... forx -> +oo 

The fitting formula (|21j) is compared with the exact BA results in Fig. |21 (top). Note that 
in the weak coupling limit our formula gives e(x — > + ) = 1/3 + x/it — 0.080x 2 , where the 
coefficient of the x 2 differs by about 4% from the exact value determined by Magyar and 
Burke by diagrammatic perturbation theory. On the other hand, the parametrization 
formula proposed by Magyar and Burke |l8| incorporates exactly this second-order weak- 
coupling term, but violates the strong-coupling asymptotic result in Eq. (|2~2|l. In fact, 
the Magyar-Burke coefficient of the 1/x term, c_i = —1.829, differs from the exact value 
c_! = -321n2/(37r) by about 22%. 

Turning to the case of attractive interactions, we find the following parametrization for- 
mula to be very accurate: 

e(x < 0) = I - ^ - A(\x\)^- (23) 
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where the function A(x), which modulates the amplitude of the strong-coupling term — x 2 /4, 
is given by 

A(x) = x ^ + a ™ x + b ™ (24) 

•E + Cjrt^C + dm 

with a m = —0.331117, b m = 0.458183, c m = a m + 4/tc, and d m = 4a m /iT + b m + 16/n 2 — 1. 
Equation (j2*3*)) embodies the exact asymptotic behaviors |24| 

. 1/3 + x/tt + ... forx^0~ , 
e(x) = { . (25) 

-x 2 /4 + 1/12 + ... forx -> -oo 

Equation ()23|) is compared with the exact BA results in Fig. El (bottom). 

The xc potential can be calculated analytically using its definition in Eq. ()12j) applied 
to Eqs. (}2~T]) or ([23)1 . The results for n GS (z) that are obtained with t£° m (n, g 1D ) determined 
according to this parametrization procedure will be termed with the acronym BALDA/2. 



III. NUMERICAL RESULTS 

We proceed to illustrate our main numerical results, which are summarized in Figs. E1E1 
In Fig. El we show the GS density profile of a Fermi gas with Af = 10 atoms at A = +2 
and —2. Repulsive interactions depress the amplitude of the shell structure, while attractive 
interactions enhance the oscillations of the profile leading to an ADW with iVf/2 distinct 
maxima related to the formation of iVf/2 spin pairs. The two BALDA schemes that we 
have proposed are in excellent agreement with each other, showing no visible difference on 
the scale of the figure. In Fig. Hjwe report a summary of our BALDA/1 results for the GS 
density profiles of a paramagnetic Fermi gas with Af = 10 atoms for increasing repulsive or 
attractive interactions. 

The shell structure is also sensitive to the system size, with larger clouds tending to have 
a relatively weaker structure. In Fig. E] we show the GS density profile of a cloud with 
Af = 20 atoms at A = +2 and —2. Comparing this figure with Fig. El it is clear that the 
amplitude of the oscillations is decreasing with increasing Af. In particular for Af = 50 
atoms at A = +2 (see Fig. El top) the Thomas- Fermi results are a good representation of 
the actual density profile, except at the edges of the cloud. On the other hand, for the same 
system size attractive interactions at A = — 2 lead to a still clearly visible ADW, though this 
is absent in the Thomas- Fermi theory (see Fig. El bottom). 

Finally, the problem of two atoms with opposite pseudospins in Q1D harmonic confine- 
ment is exactly solvable (see Appendix). A priori we do not expect an LDA approach to 
be applicable to such a small system, but from Fig. El it is seen that the BALDA scheme 
still yields some reasonable results for both repulsive and attractive interactions. For strong 
repulsive interactions it is not able to reproduce the formation of a hole at the center of the 
trap (see the discussion under point (ii) of the previous Section), while in the case of strong 
attractive interactions it overestimates the value of the density in the same region (see Fig.[7J). 
The method is nevertheless usefully reliable for the GS energy over a wide range of values of 
A (see Fig. El): for instance, at A = -30 we find £g S ALDA /(iV f fi,u;||) = -224.810 as compared to 
the exact value £ GS /(N f Hu\\) = -224.499, and at A = +30 we find £^ ALDA /(Af/?^||) = +1.01 
as compared to the exact value 8 GS /(N{hun) = +0.975. 
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IV. SUMMARY AND CONCLUSIONS 



In summary, we have presented a novel Kohn-Sham DFT study of two-component Fermi 
gases with repulsive or attractive intercomponent interactions in Q1D harmonic traps. The 
present BALDA theory, which is expected to be accurate in a weak-to-intermediate range 
of coupling strength, provides a quantitative microscopic understanding of how many-body 
exchange and correlations modify the bulk shell structure of the ground-state density profile. 
Repulsive intercomponent interactions depress the amplitude of the shell structure while 
attractive interactions stabilize atomic-density waves through Luther-Emery spin pairing. 
Such atomic-density waves should be observable in gaseous clouds containing of the order 
of up to 100 atoms, a suitable experimental technique being the search for satellites in the 
elastic diffraction pattern as discussed in Ref. [20]. It would also be important to re-examine 
numerically these GS exchange and correlation properties in relatively small systems with 
Nf < 10 atoms by exact-diagonalization or Quantum Monte Carlo methods. 

The present work can be extended in several directions. For instance, it would be in- 
teresting to generalize the present scheme to the composite-boson region B < B* in order 
to have a DFT treatment of the BCS-BEC crossover in the presence of axial confinement, 
generalizing the theories by Fuchs et al. Ill] and by Tokatly [l2|. Secondly, it would be 
interesting to study dynamical phenomena such as spin-charge separation in these strongly 
correlated gases using time-dependent DFT and/or current-DFT (§, 26|, instead of resorting 
to the inhomogeneous Tomonaga-Luttinger liquid model 2^, 23|. From a more formal DFT 



viewpoint, a functional better than in Eq. (|11|) is desirable and necessary to deal with the 
strong coupling regime. 
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APPENDIX: EXACT SOLUTION OF THE TWO- ATOM PROBLEM 

The problem of two antiparallel-spin fermions interacting with a zero-range delta-function 
potential in ID is exactly solvable thanks to the separation of centre-of-mass and relative 
variables, which is allowed by the harmonic trapping potential. Performing a canonical 
transformation to centre-of-mass [Z = [z\ + z 2 )/2,P = pi + p 2 ) and relative (z Te \ = z\ — 
Z 2,P = {pi — P2) /2) coordinates and momenta, the Hamiltonian can be written as 7i = 
Hcm(Z, P)+n ie \(z ieh p). Here, the centre-of-mass Hamiltonian H C m = P 2 / (2M)+Mu 2 Z 2 /2 
describes a free particle of mass M = 2m in a ID harmonic oscillator, while the relative- 
motion Hamiltonian 7i re i = p 2 /(2fi) + V(-2 re i) describes a free particle of mass /1 = m/2 in 
the potential V(z Te i) = ^oj 2 z 2 J2 + gm5(z re] ). 
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The spatial part of the GS wavefunction can thus be written as 

*gs(zi,z 2 ) = Mexp(-Z 2 /a})<f^(z Tel ), (A.l) 

where M is a normalization constant and (pfj(z vc i) is the GS wavefunction of the relative- 
motion problem with energy £ . Introducing the dimensionless coordinate z m \ = z Te \/a\\, the 
Schrodinger equation for the relative motion reads 



d 2 1 _ 



<rfel (*rd) = E n ^ (Zrel) , ( A.2) 



where e n = £ n /{tkO\\). Due to the antisymmetric (spin-singlet) nature of the spinorial part 
of the GS wavefunction, we need to search for the lowest (n = 0) even eigensolution of 
Eq. 

The singular delta-function term in Eq. (|A.2|) imposes a cusp on the wavefunctions at 
the origin, 

lim lim %|(M = v; ») (Sre , = 0) . (A .3) 

2 rc l-+0+ OZ re \ Zrol-^O UZ re \ 

Eq. (|A.2jl is then recognized to be the differential equation that defines the Parabolic Cylin- 
der Functions ji^]. The even solutions with the proper asymptotic behavior are |2?t l28l| 

^ ) (^rel)=^n-l/2(M), (A.4) 

Dn(x) being a Whittaker function. Using the following properties of the Whittaker func- 
tion I28j |. 



and 



together with 



^feei - 0) - 2 _ £V2+1/4r(3/4 _- n/2) (A.5) 

J™ + «^ rel 2-W2-i/4 r( i/ 4 _ fn/2 ) lA -° j 

lim = - lim (A.7) 

Zrel^O OZ Ie \ 2 rc l-»0+ OZ Ie \ 



it is easy to show that Eq. (|A.3J) becomes the trascendental equation 

r(3/4 - e n /2) A_ 

r(l/4 - e n /2) 2v/2 ' 



(A.8) 



Here r(x) is the Euler Gamma function. This equation implicitly defines the function e n (X). 
The l.h.s. of Eq. flA.8|) is shown in Fig. El as a function of — e n . The GS energy per atom is 
given by £ GB /(N{fkJ\\) = 1/4 + Eq{X)/2 and is shown in Fig. |S]as a solid line. 

Some limiting behaviors of the function Eq{\) can be established analytically from the 
properties of the Gamma function. We find 

eo(A->0) = i + -^=-£ln2 + ... (A.9) 

2 V27T 27T 
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in the weak coupling limit, 



e (A - +00) = 5 - 2 J 11 - -(ln2 - 1)- + ... (A.10) 



in the strong repulsion limit, and 



*o(A--oo) = -£-^ + A + ... (A.11) 



in the strong attraction limit. 

The GS density profile can be found from 



n GS (z)= / dz'\* GS (z,z')\ 2 . (A.12) 



The normalization constant Af is chosen according to f_ °° dzriQ$(z) = 2, i.e. 

2 3 / 2 /(v / ^ll) /A N 

N 2 = „ +QO /w v (A. 13) 

^rell^rell^rcl)! 2 



(»)/ M2_ ^>(3/4-e n /2)-^(l/4-£ n /2) 



where I2_, 

r+oc 

«,<wi 2 = — 7(1/2-17 ~ ' (A - 14) 

with -^(a;) = dlnr(x)/da;. The relation (jA.14|) holds unless e n = n + 1/2 with n = 0, 1, 2, 
when one has to use the result 

00 _ /7F 

(i2 rc i|D n (z rc i)| 2 = J-n!. (A. 15) 

For instance, in the case A = +00 we find £0 = 3/2, y^tXei) = |?rei| exp (— z 2 el /4), A/" 2 = 
2/(vra 2 ) and 

n GS W| A=+oc = [l + 2(z/a,|) 2 ] exp (-^ 2 /a 2 ) . (A.16) 

This result has form given in Eq. (|14j) . 
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FIG. 1: The dimensionless coupling parameter A as a function of the magnetic field B (in units 
of the Feshbach resonance field Bp). Here we have chosen the following values for the relevant 
parameters: m = 6.642 x 10~ 26 Kg (mass of a 40 K atom), uj± = 2ir x 100 kHz and uu = 2ir x 200 Hz 
(the anisotropy parameter of the trap is 2 x 10~ 3 ), Bp = 202.1 G, 5B = 7.8 G, and ab g = 174 Bohr 
radii. For these parameters B* = 0.991-Bf- At the 3D Feshbach resonance A F = -43.309 (see 
inset). 
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FIG. 2: Ground-state energy e GS (n, g 1D ) of the HGYM (per particle and in units of the Fermi 
energy ep) as a function of the coupling strength for a paramagnetic Fermi gas with repulsive 
interactions (top) and attractive interactions (bottom). The exact results, obtained from the 
solution of the BA equations are compared with the fitting formulae in Eq. (|21j) and in 

Eq. (solid lines). 



17 




FIG. 3: Density profile ugs( z ) (i n units of ) as a function of z/a\\ for a paramagnetic Fermi gas 
with Nf = 10 atoms at A = +2 and —2. The results of the BALDA/1 scheme are compared with 
those of the BALDA/2 scheme. The thin solid line corresponds to the noninteracting A = case. 
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FIG. 4: Evolution of the density profile Uqs(z) (in units of o^ 1 ) with increasing A in the BALDA/1 
scheme, for a paramagnetic Fermi gas of iVf = 10 atoms with repulsive interactions (top) and 
attractive interactions (bottom). The curve at A = +oo is the theoretical result given in Eq. (|14[). 
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FIG. 5: Density profile riQs(z) (in units of Q^ 1 ) a s a function of z/a^ for a paramagnetic Fermi 
gas with Nf = 20 atoms at A = +2 and —2. The thin solid line corresponds to the noninteracting 
A = case. 
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FIG. 6: Density profile tiqs(z) (in units of Oil 1 ) as a function of z/a\\ for paramagnetic Fermi gases 
with JVf = 20 and 50 atoms at A = +2 (top) and A = -2 (bottom). The results of the BALDA/1 
scheme are compared with the Thomas- Fermi results. 
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FIG. 7: Density profile nQs(z) (in units of fly" 1 ) as a function of z/a^ for two Fermi atoms with 
opposite pseudospins at A = +1 and — 1 (top), A = +2 and —2 (middle), and A = +10 and —10 
(bottom). The results of the BALDA/1 scheme are compared with the exact results. 
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FIG. 8: Ground-state energy per atom E GS /N{ (in units of TkJn) as a function of A for two Fermi 
atoms with opposite pseudospins. The results of the BALDA/1 scheme are compared with the 
exact results. 
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FIG. 9: A plot of the function h{x) = r(3/4 + x/2)/r(l/4 + x/2). h(x) has zeroes at x = -2n-l/2 
and poles at x = — 2n — 3/2, with n = 0, 1,2, .... In order to find the GS energy of the 2-atoms 
system one has to find the intersections of horizontal lines with the branch of h(x) in the range 
—3/2 < x < +oo. The upper (lower) half-plane is relevant for attractive (repulsive) interactions. 
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